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There has been rapid development of systems that yield strong interactions between freely propa¬ 
gating photons in one dimension via controlled coupling to quantum emitters. This raises interesting 
possibilities such as quantum information processing with photons or quantum many-body states of 
light, but treating such systems generally remains a difficult task theoretically. Here, we describe a 
novel technique in which the dynamics and correlations of a few photons can be exactly calculated, 
based upon knowledge of the initial photonic state and the solution of the reduced effective dynamics 
of the quantum emitters alone. We show that this generalized “input-output” formalism allows for 
a straightforward numerical implementation regardless of system details, such as emitter positions, 
external driving, and level structure. As a specific example, we apply our technique to show how 
atomic systems with infinite-range interactions and under conditions of electromagnetically induced 
transparency enable the selective transmission of correlated multi-photon states. 


I. INTRODUCTION 

Systems in which individual photons can interact 
strongly with each other constitute an exciting frontier 
for the fields of quantum and nonlinear optics [lj. Such 
systems enable the generation and manipulation of non- 
classical light, which are crucial ingredients for quantum 
information processing and quantum networks [2]. At 
the many-body level, it has been predicted that these 
systems can produce phenomena such as quantum phase 
transitions of light om or photon crystallization 00. 
Early examples of such systems where strong interactions 
between photons could be observed consisted of individ¬ 
ual atoms coupled to single modes of high-finesse opti¬ 
cal cavities, within the context of cavity quantum elec¬ 
trodynamics (QED) [7, 8]. More recently, a number of 
systems have emerged that produce strong optical non- 
linearities between freely propagating photons, including 
cold atomic gases coupled to guided modes of tapered 
fibers mm, photonic crystal fibers El or waveguides 
[12] , cold Rydberg gases in free space [13j , and supercon¬ 
ducting qubits coupled to microwave waveguides mm- 

The paradigm of cavity QED admits exact analytical 
or numerical solutions at the few-photon (and/or few- 
atom) level through the elegant “input-output” formal¬ 
ism m, as illustrated conceptually in Fig. [lja). In this 
formalism, all of the properties of the field exiting the sys¬ 
tem (the output) can be determined based upon knowl¬ 
edge of the input field and the dynamics of the atom- 
cavity system alone. For a few excitations, the latter 
can be solved due to the small Hilbert space associated 
with a small number of excitations of a discrete cavity 
mode. In contrast, despite rapid development on the 
experimental front, theoretical techniques to treat the 
dynamics of freely propagating photons interacting with 
spatially distributed emitters are generally lacking. At 
first glance, the challenge compared to the cavity case 
arises from the fact that a two-level system is a non¬ 
linear frequency mixer, which is capable of generating a 


continuum of new frequencies from an initial pulse, as 
schematically depicted in Fig. 0b). A priori, keeping 
track of this continuum as it propagates and re-scatters 
from other emitters appears to be a difficult task. An ex¬ 
ception is the weak excitation limit, in which atoms can 
be treated as linear scatterers and the powerful transfer 
matrix method of linear optics can be employed Bang. 

The full quantum case has been solved exactly in a lim¬ 
ited number of situations in which nonlinear systems are 
coupled to ID waveguides, such as up to three photons 
scattering on a single atom mm, one and two pho¬ 
tons scattering on a cavity QED system m, and many 
photons scattering on many atoms coupled to a chiral 
(monodirectional) waveguide [22 . The formalism em¬ 
ployed in Ref. [2d is particularly elegant, because it es¬ 
tablishes an input-output relation to determine the non¬ 
linear scattering from a two-level atom. Here, we show 
that this technique can be efficiently generalized to many 
atoms, chiral or bi-directional waveguides, and arbitrary 
atomic configurations, providing a powerful tool to inves¬ 
tigate nonlinear optical dynamics in all systems of inter¬ 
est. 

This paper is organized in the following way: first, 
we present a generalized input-output formalism to treat 
few-photon propagation in waveguides coupled to many 
atoms. We show that the infinite degrees of freedom 
associated with the photonic modes can be effectively in¬ 
tegrated out, yielding an open, interacting “spin” model 
that involves only the internal degrees of freedom of the 
atoms. This open system can be solved using a number 
of conventional, quantum optical techniques. Then, we 
show that the solution of the spin problem can be used to 
re-construct the optical fields. In particular, we provide 
a prescription to map spin correlations to S-matrix ele¬ 
ments, which contain full information about the photon 
dynamics, and give explicit closed-form expressions for 
the one- and two-photon cases. Importantly, in analogy 
with the cavity QED case, our technique enables analyt¬ 
ical solutions under some scenarios, but in general allows 
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for simple numerical implementation under a wide variety 
of circumstances of interest, such as different level struc¬ 
tures, external driving, atomic positions, atomic motion, 
etc. Finally, to illustrate the ease of usage, we apply 
our technique to the study of nonlinear field propaga¬ 
tion through an optically dense ensemble of atoms with 
Rydberg-like interactions. 

II. GENERALIZED INPUT-OUTPUT 
FORMALISM 

In this section we consider a generic system composed 
of many atoms located at positions zi along a bidirec¬ 
tional waveguide. We assume that there is an optical 
transition between ground and excited-state levels | g) 
and \e) to which the waveguide couples, but otherwise 
we leave unspecified the atomic internal structure and the 
possible interactions between them ( e.g ., Rydberg inter¬ 
actions), as such terms do not affect the derivation pre¬ 
sented here. The bare Hamiltonian of the system is com¬ 
posed of a term describing the energy levels of the atoms 
H at , and a waveguide part H ph = Y. v =± f dkio k b\ k b v ^, 
where k is the wavevector and v = =b is an index for the 
direction of propagation, with the plus (minus) denot¬ 
ing propagation towards the right (left) direction. We 
assume that within the bandwidth of modes to which 
the atoms significantly couple, the dispersion relation for 
the guided modes can be linearized as ujk = c|fc|. The 
interaction between atoms and photons is given by 

N 

tfint=5EE 

i= 1 
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FIG. 1 . Cavity QED vs. many-atom waveguide (a) 

Schematic representation of a cavity QED system. The input- 
output formalism enables the outgoing field to be calculated 
based on knowledge of the input field and the internal dy¬ 
namics of the cavity system, (b) System composed of many 
atoms coupled to a common waveguide. The nonlinearity of 
an atom enables the generation of a continuum of new fre¬ 
quencies upon scattering of an incoming state. This property, 
combined with multiple scattering from other atoms, appears 
to make this system more complicated than the cavity QED 
case. 

of the atomic transition can be retained by introducing 
an interaction energy Uo for multiple excitations, through 
the term ( Uq/2 ) aja*(aja* — 1) in i7 at , and by taking 
the limit Uo oo at the end of calculations for observ¬ 
able quantities. 

The Heisenberg equations for b v ^ can be formally inte¬ 
grated and Fourier transformed, to obtain the real-space 
wave equation 


/ 


dk (b Vtk a l eg e lvkZi + h.c.), (1) 


which describes the process where excited atoms can emit 
photons into the waveguide, or ground-state atoms can 
become excited by absorbing a photon. The coupling 
amplitude g is assumed to be identical for all atoms, 
while the coupling phase depends on the atomic position 
( ( elvkZiy Here, we will explicitly treat the more compli¬ 
cated bidirectional case, although all of the results readily 
generalize to the case of a single direction of propagation. 

Our immediate goal is to generalize the well-known 
input-output formalism of cavity QED m and the more 
recent formalism for a single atom coupled to a waveguide 
m , to the present situation of many spatially distributed 
atoms coupled to a common waveguide. In short, we will 
eliminate the photonic degrees of freedom by formal in¬ 
tegration, which reveals that the output field exiting the 
collection of atoms is completely describable in terms of 
the input field and atomic properties alone. This formal 
integration also provides a set of generalized Heisenberg- 
Langevin equations that governs the atomic evolution. 

The Heisenberg equations of motion for a ge and b v ^ 
can be readily obtained by calculating the commutators 
with H. To simplify the presentation of the resulting 
equations, we replace the spin operators a ge with the 
bosonic annihilation operator cq. The two-level nature 


6 v (z,t) — hy 5 in(t vz/c) 



E 8{v(z - Zj))dj (t - V ( z - Zj U p) 


Here b v j n is the homogeneous solution, physically corre¬ 
sponding to the freely propagating field in the waveguide, 
while the second term on the right consists of the part 
of the field emitted by the atoms. Inserting Eq. © into 
the equation for cq, we obtain 


di = i[H at , di] - i V2ng E b v,m(t - vzt/c) 




27 rg 2 


N 


E aj(t- \zi~Zj\/c)). 


(3) 


3 = 1 


In realistic systems, time retardation can be neglected, 
resulting in the Markov approximation aj(t — \zi — 
Zj |/c) ~ aj{t)e luJ '^ Zi ~ Zj ^ c . Here, u;i n is a central fre¬ 
quency around which the atomic dynamics is centered 
(typically the atomic resonance frequency). This approx¬ 
imation is valid when the difference in free-space propa¬ 
gation phases Ac oL/c <C 1 is small across the characteris¬ 
tic system size L and over the bandwidth of photons Ac o 
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involved in the dynamics. As a simple example, the char¬ 
acteristic bandwidth of an atomic system is given by its 
spontaneous emission rate, corresponding to a few MHz, 
which results in a significant free-space phase difference 
only over lengths L > 1 m much longer than realistic 
atomic ensembles. 

We have thus obtained the generalized Heisenberg- 
Langevin equation 


hi — ^\Z^rip/2 ^ ^ b v \ n (t VZi/c ) 


v=± 


v N 


(4) 


3 = 1 


where we have identified T id = 47rg 2 /c as the single¬ 
atom spontaneous emission rate into the waveguide 
modes. If we keep separated the terms proportional to aj 
coming from the right and left-going photonic fields, we 
can find easily that the Lindblad jump operators corre- 
spondin g to the decay of the atoms into the waveguide are 
0± = ^/Fid/ 4 (ije TluJiIlZ i/ c : in terms of which we can 
write the master equation for the atomic density matrix 

P = £{P\ = ,P\+Y,V = ± 2 0 V pOl~OlO V p-pOlO V . 

We also see that we can derive Eq. @ from a non- 
Hermitian effective Hamiltonian 

H eS = H at - ^ Z e^~ z ^ c , (5) 

ij 


which can be used for a quantum jump description of 
the atomic dynamics. The resulting infinite-range in¬ 
teraction between a pair of atoms ij intuitively results 
from the propagation of a mediating photon between that 
pair, with a phase factor proportional to the separation 
distance. 

Within the same approximations employed 
above to derive the Heisenberg-Langevin equations 
we can obtain a generalized input-output rela¬ 
tion of the form b v , out (z,t) = b v ^ n (t — vz/c) — 
iy/T id /(2c) ai(t) e lvuJin ( z ~ Zi ^ c , where the output 

field is defined for z > zr = max[^] (z < zr = min [ 2 ^]) 
for right (left)-going fields. However, since the right-going 
output field propagates freely after zr , it is convenient 
to simply define fr+, 0 ut(£) = fr+,out(^R + e,£) as the field 
immediately past the right-most atom (where e is an 
infinitesimal positive number), and similarly for the 
left-going output. The derived relation shows that the 
out-going field properties are obtainable from those of 
the atoms alone. 

The emergence of infinite-range interactions between 
emitters mediated by guided photons, and input-output 
relationships between these emitters and the outgoing 
field, have been discussed before in a number of con¬ 
texts muni ns, but the idea that such concepts could 
be used to study quantum interactions of photons in ex¬ 
tended systems has not been fully appreciated. In the 
remaining sections, we will demonstrate the effectiveness 


of this approach to quantum nonlinear optics. In partic¬ 
ular, the infinite-dimensional continuum of the photons 
is effectively reduced to a Hilbert space of dimension 
dim^] = Ci) where n is the maximum number 

of atomic excitations (for n = TV we have dim [H\ = 2^). 
The atomic dynamics, on the other hand, having been re¬ 
duced to standard Heisenberg-Langevin equations, quan¬ 
tum jump, or master equations, are solvable by conven¬ 
tional prescriptions [24] . It is also possible to derive gen¬ 
eralized master equations describing the atomic dynam¬ 
ics in response to arbitrary (e.g., non-classical) incident 
states of light, as detailed in Appendix C. 


III. RELATION TO S-MATRIX ELEMENTS 

The S-matrix characterizes how an incoming state 
of monochromatic photons evolves via interaction with 
atoms into a superposition of outgoing monochromatic 
photons. Because monochromatic photons form a com¬ 
plete basis, the S-matrix thus contains all information 
about photon dynamics. Here, we will show that S- 
matrix elements can be calculated from the dynamics of 
the atoms evolving under the effective spin model. 

Formally, the S-matrix for the interaction of an 71 - 
photon state with an arbitrary system is defined as 

S { p n l = (p\S\k) 

= (0| ^out (pi) -•bout (p n )bUki)A a (k n ) |0> 

= <0| b out (h).A ut (t n )bi(t[)..bl n (0 |0)(6) 

where the input and output creation operators respec¬ 
tively create freely propagating incoming and outgo¬ 
ing photonic states. The vectors p and k denote 
the outgoing and incoming frequencies of the n pho¬ 
tons. The input and output operators can be any 
combination of + and - propagation directions (we 
have omitted this index here for simplicity). In the 
last line we have used a global Fourier transformation 
jrj(2n) = (27r)- n /nr=i cM*ie ic(tiP4_t < fe< ), to express 
the S-matrix in time. 

The S-matrix has been calculated exactly before in 
a limited number of situations fT9ff22] . Here, we pro¬ 
vide a general prescription to numerically obtain the S- 
matrix starting from the input-output formalism. First, 
it should be noted that the operators 6i n , b Q ut correspond 
to the input and output operators defined in the previous 
section m ■ On the other hand, the input-output rela¬ 
tion enables the correlator of Eq. © to be written purely 
in terms of atomic operators. For notational simplicity, 
we give a derivation for a single spin and a monodirec- 
tional waveguide, but its generalization to the bidirec¬ 
tional waveguide and many atoms is straightforward. For 
our purpose it is enough to have an input-output relation 
of the form 


front = frin - i\/7 a : 


(7) 
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where a is in our case the spin operator a ge . 

We summarize the main idea of the derivation here, 
while the details can be found in Appendix A.l. We be¬ 
gin by noting that Eq. 0 enables one to replace output 
operators by a combination of system and input opera¬ 
tors, or input operators by system and output operators. 
Selectively using these substitutions, one can exploit fa¬ 
vorable properties of either the input or output field, in 
order to gradually time order all of the system opera¬ 
tors (where operators at later times appear to the left of 
those at earlier times), while removing input and output 
operators from the correlation. The favorable properties 
that can be used are that 1) system and input operators 
commute, [a(t),b{ n (t f )] = 0, when t' > £, and likewise 
[a(t), b out (t')] = 0 for t' < £, 2) input annihilation opera¬ 
tors at different times commute amongst themselves, as 
do output annihilation operators, and 3) the correlator 
can be reduced in size using [bi n (t),bl n (t')] = S(t — £'), 
or made to vanish using 6 in (t) |0) = 0 or (0| bl ut (t) = 0. 
Through this procedure, the S-matrix can be expressed 
as a Fourier transform of a sum of terms involving only 
time-ordered atomic operators (indicated by the operator 
T). These functions generally have the form 

(0\T\a(t 1 )..a(t m )a\t' 1 )..a\t' m )] |0), (8) 

with m < n, multiplied by n — m delta functions in time. 

Moreover, using the general expression for the 
Heisenberg-Langevin equation and the quantum regres¬ 
sion theorem, it can be proven that when external fields 
driving the system do not generate waveguide photons, 
the correlation function of Eq. (|8| can be evaluated by 
evolving a(t) as e lHef{t a e~ lHe{ft (see Appendix A.2 for a 
formal derivation). Here, H e ff is the effective Hamilto¬ 
nian from Eq. 0 that contains only the spin operators. 
Although such a form for a(t ) is not true in general due 
to quantum noise, these noise terms have no influence 
on the correlation. A similar procedure as above enables 
one to express other important observables of the field, 
such as the second-order correlation function g^ 2 \t), in 
terms of correlation functions involving atoms alone. 

While the discussion has thus far been completely gen¬ 
eral, the case of S-matrix elements involving only one 
or two photons can be formally reduced to particularly 
simple expressions. For example, in Appendix B, we 
show that the transmission coefficient X \ for the many- 
atom, bi-directional waveguide case is related to the S- 
matrix by S^l. k+ = (0| b +tOUt (p)b\_ in (k) |0) = T k 5 pk . 
Furthermore, it can be expressed in terms of a known 
~ N x N matrix corresponding to the single-excitation 
Green’s function Go (whose form varies depending on the 
system details), 

T k = 1 - ^ E [<*»(*)]« (9) 

ij 

Similarly, the two-photon S-matrix in transmission is 





FIG. 2. (a) EIT level scheme. The atomic ground (|g)) and 
excited states (|e)) interact with the quantum propagating 
field b of the waveguide. An additional classical field with 
Rabi frequency Ft couples state |e) to a metastable state | s). 
The total single-atom linewidth of the excited state is given 
by r. (b) The real (y 7 ) an d imaginary (y") parts of the linear 
susceptibility for a two-level atom (upper panel) and three- 
level atom (lower panel), as a function of the dimensionless 
detuning S/T of the field b from the resonance frequency of the 
\g)-\e) transition. For the three-level atom, the parameters 
used are Sl = 0 and Fl/F = 1/3. 


generally given by 


o( 2 ) _ 

°Pl + ,P2 + ’M+,k 2 + 

F 2 

_ \ 1 A 

1 u Pi+P2M+k 2 


Tk 1 Tk 2 Spik i 3p 2 k 2 T 

WijTij-i'j’Wi'j' + (pi <-> P2), 

iji'j' 


( 10 ) 


where the first and second terms on the right describe the 
linear and nonlinear contributions, respectively. The lat¬ 
ter term can be expressed in terms of matrices W related 
to the single-excitation Green’s function, and a known 
~ N 2 x N 2 matrix T characterizing atomic nonlineari- 
ties and interactions. 


IV. ELECTROMAGNETIC INDUCED 
TRANSPARENCY 

In this section, we apply our formalism to a spe¬ 
cific example involving three-level atoms under condi¬ 
tions of electromagnetically induced transparency (EIT) 
and with Rydberg-like interactions between atoms m- 
The linear susceptibility for a two-level atom with states 
| g) and |e), in response to a weak probe field with de¬ 
tuning 6 = cOp — uj eg from the atomic resonance, is shown 
in Fig. It- It can be seen that the response on reso¬ 
nance is primarily absorptive, as characterized by the 
imaginary part of the susceptibility (y", red curve). In 
contrast, the response can become primarily dispersive 
near resonance if a third level | s) is added, and if the 
transition \e) — \s) is driven by a control field (charac¬ 
terized by Rabi frequency Ft and single photon detuning 
Sl = — cj es ). Specifically, via interference between 












5 



FIG. 3. EIT: single polariton propagation for a p < <jeit (a) 
and a p > <jeit (b). Plotted is the population Pj = ( a J ss ) 
of atom j in the state |s). The blue line corresponds to the 
initial state, the green line to the state evolved over a time 
tf , and the red dots to the theoretically predicted evolution. 
Other parameters: tfV g = Nd/6, T = 5, Y\d — 10, Q = 1, 
N = 500 and a eit ~ 22 d, where d is the lattice constant. 

the probe and control fields, the medium can become 
transparent to the probe field (x" = 0) when two-photon 
resonance is achieved, S — Sl =0, realizing EIT [25] . 
In this process, the incoming probe field strongly mixes 
with spin wave excitations a sg to create “dark-state po- 
laritons.” The medium remains highly transparent within 
a characteristic bandwidth A eit around the two-photon 
resonance, which reduces to A eit ~ 2 Q 2 /(Ty/~D), when 
Sl = 0. Here we have introduced the total single-atom 
linewidth T (see below) and the optical depth D , cor¬ 
responding to the opacity of the medium in absence of 
EIT and defined in terms of input and output intensity 
as / out (fi = 0)//in(£2 = 0) = exp(— D). These polaritons 
propagate at a strongly reduced group velocity v g <C c, 
as indicated by the steep slope of the real part of the 
susceptibility x' in Fig. [ 2 |), which is proportional to the 
control field intensity [25] . 

Taking Si to be the annihilation operator for the state 
| Si), EIT is described within our spin model by the effec¬ 
tive spin Hamiltonian 

p/ _ _ 

H = —(Sl + 1—) Y a ) a i ~ Q + s i a i) 

3 3 

-N^Y eikinlZj ~ zlla N (n) 

3,1 

where the first line represents the explicit form of H at in 
Eq. © for the EIT three-level atomic structure. In ad¬ 
dition to waveguide coupling, here we have added an in¬ 
dependent atomic decay rate T' into other channels (e.g., 
unguided modes), yielding a total single-atom linewidth 

ofr = r + r 1D . 

Our theoretical model nearly ideally describes exper¬ 
iments in which atoms are coupled to one-dimensional 
waveguides such as nanofibers [9] or photonic crys¬ 
tals [TO], in which the interaction probability of a single 
photon and atom is T m/T ~ 0.1 is significant and up 
to A - 10 3 atoms are trapped. One consequence of the 


large interaction probability is that even a single atom 
can yield a significant reflectance for single photons m • 
In free-space experiments, the atom-photon interaction 
probability is much smaller, while a much larger number 
of atoms is used to induce a significant optical response. 
While this large number cannot be implemented numeri¬ 
cally with our model, our system can nonetheless be used 
to reproduce macroscopic observables, thus making our 
approach an excellent description for atom-light inter¬ 
faces in general. 

Intuitively (and as can be rigorously shown, see below), 
provided that free-space experiments only involve a single 
transverse mode of light, one expects that the system re¬ 
sponse to light only depends on the interaction probabil¬ 
ity and number of atoms via their product. This product 
in fact defines the optical depth D = N(2T 1D /T) [26] . 
As D < 10 2 in free-space experiments, such systems can 
be evaluated using our spin model with several hundred 
atoms and suitably large Tid/T. The only remaining is¬ 
sue is the potentially large reflection generated by atoms 
in the mathematical model, compared to free-space en¬ 
sembles where reflection is negligible. Reflection can be 
suppressed in the model by giving the atoms a lattice 
constant d where k in d = (2 m + l)7r/2, where m is a non¬ 
negative integer, such that the reflection from atoms de¬ 
structively interferes [29, 30] . All subsequent numerical 
calculations are performed under this condition. 

The spin model on a lattice describing EIT, Eq. 
can be exactly solved in the linear regime using the trans¬ 
fer matrix formalism [29] , which correctly reproduces the 
free-space result and dependence on optical depth for 
group velocity v g = 2f2 2 n/Ti£> and transparency win¬ 
dow A eit, where n is the (linear) atomic density. The 
corresponding minimum spatial extent of a pulse that can 
propagate inside the medium with high transparency is 
given by cr eit = v g / Aeit - 

A single photon propagating inside an ensemble of 
atoms under EIT conditions is coherently mapped onto 
a single dark polariton, corresponding to a delocalized 
spin wave populating the single excitation subspace of 
the atomic ensemble. The polariton dynamics can be 
therefore visualized directly by monitoring the excitation 
probability (c t j ss ) of the atoms in the ensemble. In Fig. [ 3 J 
we initialize a single polariton inside the medium with an 
atomic wave function of the form | x/j) = JA fj& J sg \g)® N , 
and we determine numerically the time evolution under 
H in Eq. ( [IT] ) up to a final time tf. Choosing an ini¬ 
tially Gaussian spin wave, fj = exp(i/q n <ijf)exp(—(j<i — 
/i) 2 /4<j 2 )/(27T(j 2 ) 1 / 4 , with spatial extent cr p (blue line), 
one sees that the wavepacket propagates a distance v g -tf , 
and with little loss provided that a p > cfeit- Numer¬ 
ics (green line) show perfect agreement with theoretical 
predictions (red lines) obtained via the transfer matrix 
formalism m ■ 
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FIG. 4. (a) single-photon (circles) and two-photon (trian¬ 

gles) transmission spectrum for a weak probe field, for se¬ 
lected values (7/2 = 0 and (7/2 = 0.2 of the infinite-range 
interaction strength. The linear transmission is independent 
of (7, while the two-photon spectrum exhibits a shift in the 
maximum transmission by an amount (7/2. Other parame- 
ters: N = 200, r 1D = 1,T' = 3, « = 2, Sl = 0, £ = 10“ 6 . 
(b) Contour plot of the two-photon transmission spectrum 
T 2 = (b\ 

,out (t)b\ 

,out (t)b+ ,out (t)b+ ,out (■ t))/S 4 , as functions of 
interaction strength (7/2 and two photon detuning A — S 
(Sl = 0). Cuts of the contour plot (illustrated by the dashed 
lines) are plotted in a), (c) Comparison of g^ 2 \r) evalu¬ 

ated by numerical simulations (red dashed line) and S-matrix 
theory (black line). For constant infinite range interactions 
(7=1 and A = 0, the interactions induce photon antibunch¬ 
ing. Other parameters: N = 20, Fid = 2 , T' = 2 , Q = 1, 
S L = 0, S = 1CT 6 . 


V. INFINITE RANGE INTERACTION 

The spin model formalism can be easily extended to 
include arbitrary atomic interactions, providing a power¬ 
ful tool to study quantum nonlinear optical effects. As a 
concrete example, we consider a system in which atoms 
can interact directly over a long range, such as via Ryd¬ 
berg states or photonic crystal bandgaps [33] . 

The total Hamiltonian is given by 

p/ _ _ 

H = ~{8 l + i—) ^2 a ) a i ~ n $7 a i s i + s ] a i ) 

3 3 

~ E e ik ^-^a] ai + \Y1 Ui jS }sjs lSj 

3,1 3,1 

+ ^drive (12) 

in which Uij represents a dispersive interaction between 
atoms i and j when they are simultaneously in state 
| s). As we are primarily interested in demonstrating the 
use of our technique, we take here a “toy model” where 
atoms experience a constant infinite-range interaction, 
Uij /2 = (7. Such a case enables the numerical results to 
be intuitively understood, although we note that other 
choices of Uij do not increase the numerical complexity. 
In particular we are interested in studying the propa¬ 
gation of a constant weak coherent input field through 
the atomic ensemble. The corresponding driving then is 
given by i^drive = £ ^2j(aje lkinZj e~ lAt + aje~ lkinZj e lAt ), 
where £ « T is the amplitude of the constant driving 
field, A = S — Sl the detuning from two photon res¬ 
onance condition, and the initial state is given by the 
global atomic ground state |^) = \g)® N [27j |28] . With 
infinite-range interaction, one spin flip to state | s)j shifts 
the energies of all other states | s)j by an amount (7. A 
second photon should then be able to propagate with per¬ 
fect transparency, provided it has a detuning compensat¬ 
ing for the energy shift (7, thus ensuring the two-photon 
resonance condition is satisfied. As we result, we expect 
to see a transparency window for two photons, whose 
central frequency shifts linearly with (7. 

This predicted behavior can be confirmed 
by plotting the transmitted intensity frac- 
tion, Ti = //Tin = (<&+,out (*)&+,°ut (*))/£ 2 , and also 
the second-order correlation function T\ = 
( & +,out( i ) & +,out( i ) & +,out(i)&+,out(i))/f 4 , Which corre- 
sponds roughly to the two-photon transmission. Fig. [4] 
shows the single-photon transmission T\ and two-photon 
transmission T\ as a function of the interaction strength 
(7 and detuning from two photon resonance A. As 
expected, T\ shows a peak at A = 0 independently of 
the interaction intensity (7; instead the peak in T<± shifts 
towards A = (7/2 with increasing (7. The decay of X \ 
for increasing (7 can be intuitively understood by noting 
that we have a constant coherent state input, in which 
photons are randomly spaced, causing two photons to 
enter the medium at different times. Thus, until the 
second photon enters, the first photon propagates as a 
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single polariton detuned by A from the single-photon 
transparency condition, getting partially absorbed in 
the process. By increasing the interaction we increase 
the detuning for this single polariton and consequently 
its absorption, explaining the trend observed for T 2 in 
Fig.@ A quantitative description of this phenomenon is 
given in Appendix D. 

Field correlation functions like intensity 

1 = A,ouiV)kf,out (t)) or g 2 (r) = (&touiV) & +,ouiV + 
r)b+, out (t + r)b +) out 0 m 2 can be computed ac¬ 
cording to the following strategy. First we switch 
from Heisenberg representation to Schrodinger 
representation, so that for the intensity we get: 

1 = ( & +,out(*)&+,out(*)) = (V’(i)l&+,out & +,out|V’(i))- The 
time evolved wave function, |is determined by nu¬ 
merically evolving the initial spin state |^) under H for 
a time t. Then, the state immediately after detection of 
one photon, & +jOUt |'0(t)} = 10), is evaluated by expressing 
b +5out in terms of spin operators using the input-output 
formalism: b +5ou t = £e lkinZR — iriz>/2 JA a ge e lkin ^ ZR ~ Zj \ 
Finally we obtain the intensity by computing the 
probability of the one-photon detected state, I = (0|0). 


For g 2 (t) an extra step is needed. Its numerator 
describes the process of detecting two photons, the 
first at time t and the second at time t + r: f(t + 

t ) = (&+,outWfr+,out(£ + r)b +tOU t(t + r)b +)OUt (i)). As 
before we can switch to Schrodinger picture, f(t- K 

r) = ('4’(t)\b\_ out e iHT 6 + ;OU t&+,oute -iffT b +tOUt \ip(t)), and 
evaluate the state after dectection of the first photon, 
fr+,outlVKO) = 10)- Then detection of a second pho¬ 
ton after a time r entails performing an extra evo¬ 
lution under H and annihilating a photon, that is: 
6+,oute -1 ^ r |0) = b+,out|0( T ))- Finally, we evaluate the 

quantity f{t + r) = (0(r)|&t, O ut & +,out|<K'r))> by again ex- 
pressing b +?out in terms of spin operators. In Fig. [4^, 
we plot the numerically obtained result for g( 2 \r), for 
the case where infinite-range interactions are turned on 
((7 = 1) and for a weak coherent input state with detun¬ 
ing A = 0. In such a situation, one expects for the single¬ 
photon component of the coherent state to transmit per¬ 
fectly, while the two-photon component is detuned from 
its transparency window and becomes absorbed. This 
nonlinear absorption intuitively yields the strong anti¬ 
bunching dip g^ 2 \r = 0) < 1. We also evaluate this 
second-order correlation function using the analytical re¬ 
sult for the two-photon S-matrix in Eq. (10), which shows 
perfect agreement as expected. 


VI. CONCLUSION 

We have shown that the dynamics of a few photons, 
propagating under strong interactions mediated by quan¬ 
tum emitters, is fully and efficiently characterized by the 
dynamics of an open quantum spin model. As the spin 
model is solvable by standard quantum optical techniques 


for open systems, our approach provides an easily im- 
plementable recipe for the exact numerical study of a 
large class of quantum nonlinear optical systems. As 
an example, it provides an attractive alternative to nu¬ 
merical simulations of propagation through cold Rydberg 
gases, where typically the continuous field is finely dis¬ 
cretized and its degrees of freedom are explicitly kept 
track of m- Our technique can also be used to treat a 
number of strongly nonlinear systems based upon atomic 
saturation |[5j, where only approximate effective theories 
for field propagation were previously available. We an¬ 
ticipate that the availability of exact results will greatly 
aid in the development of effective theories for the gener¬ 
ally challenging problem of quantum many-body states of 
light, and that our model will shed new insight on con¬ 
ceptual links between such systems and quantum spin 
systems. 
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Appendix A: Scattering theory 

1. Connecting S-matrix to atomic correlation 
functions 

We present here the full derivation of the decomposi¬ 
tion of the S-matrix elements in terms of time-ordered 
atomic correlation functions. The starting point is the 
definition of the S-matrix in Eq. ([6|. Since the output 
operators commute between themselves because of the in- 
distinguishability of photons, they can be freely ordered 
by decreasing times. Introducing the time ordering op¬ 
erator T and also using Eq. (|7|), the operators in Eq. (J6| 
can be written as 

T[(b in (ti) + y/^a(ti))..(b in (t n ) + V7 a{t n ))] x 

x 4Vl)4n(0- ( A1 ) 

It is natural to label the terms above by the number 
m of system operators a present in each term. Thanks 
to the fact that [a(t), bin(O] = 0 for t' > £, all the input 
operators can be moved to the right of the spin operators. 
Thus, the term of order m will be of the form 

(0\T[a(t 1 )..a(t m )] bi n (t m+ i)..b in (t n )bl n (t[)..bl n (t' n ) |0), 

(A2) 

and can be simplified using the commutation relations 
between input operators [bi n (t),bl n (t')} = S(t — t'). This 



manipulation results in a sum of (n!) 2 /(m!) 2 (n — m)\ 
terms for each original term of order m. Each term of 
the sum consists of n — m delta functions multiplied by 
a correlation function of the form 

(0|T[o(t 1 )..a(t m )]6j n (ti)..6| I 1 (C) |0>. (A3) 

Since the a operators commute with the b\ n operators at 
later times, the time ordering operator can be extended 
to all the operators in the correlation function. Using 
again Eq. ([T]) to express the input operators one gets 

(0|T[a(ti)..a(t m ) - ^ 7 a\t[))..x 

x (&L(C)-V7« f (C))]|0). (A4) 

However, since the operators bl ut commute with all the 
operators a on the left, only 

(0| T[a(t 1 )..a(t m )a\t' 1 )..a\t' m )\ |0) (A5) 

remains, which reproduces Eq. <©• 

2. Evolution under the effective Hamiltonian 


1. n photons scattering on a single two-level atom 

The S-matrix element for the scattering of n pho¬ 
tons can be in general decomposed in a part that is a 
series of products of lower-order elements and a part that 
cannot be expressed in such a way. The latter is called 
fully connected part of the matrix element and physically 
corresponds to a n-body interaction, and is denoted by 
iT( n ). Thus, knowing how to decompose the S-matrix, 
one has to calculate iwith 1 < j < n to construct 
S^ n \ Furthermore, it can shown that the fully connected 
part in an element of order n can be obtained by the sys¬ 
tem correlation function of order n. We show here how to 
calculate such elements using the results presented above, 
for the case of a single two-level atom coupled to a one- 
directional waveguide. It will be possible to appreciate 
the simplicity of our formalism compared the more cum¬ 
bersome method used in Ref. [ 22 ] to obtain the same 
result. 

We start from the relation between the connected part 
of the matrix element and the correlation function of 
atomic operators 


In Eq. (J 8 | the vacuum state |0) stands for |0 ) 6 \g)® N , 
i.e., the vacuum state of field modes and the ground state 
of all the atoms. The atomic operators are in the Heisen¬ 
berg picture, a(t) = e lHt ae ~ lHt , and H is the Hamilto¬ 
nian of the whole system without the driving field. By 
taking the term with t\ > £ 2 ... > t m > > t f m as 

an example (our argument holds for any time ordering), 
we show now that Eq. can be evaluated by effectively 
evolving system operators as a(t) = e ljHefft ae _ljHefft . The 
quantum regression theorem is applied here to eliminate 
the bath or photonic degree of freedom, and results in 

< 0| 6 ( 0 | a(t 1 )...a(t m )a' ! (t' 1 )...a\t' m ) | 0) 6 | 0 ) 

= Tr[ae c( - tl ~ t2 ' ) a...a* e c ( tm - 1 ~ tm ' ) at p( 0 )], (A 6 ) 



x (Ta-(t n )...a-(t 1 )a+(t' n )...a + (t , 1 )), (Bl) 

where a _(i) = e lHettt cr—e~ lHettt , with H e g = (oj eg — 
iT/2)cr ee . The time ordering in the correlator gives 2n! 
possible orderings, but it is easy to arrive to the conclu¬ 
sion that only orderings which start on the left with a cr_ 
operator and alternate cr + and cr_ give a non-zero con¬ 
tribution. It is also immediate to see that the number of 
this possible orderings is (n!) 2 . A possible ordering is for 
instance 

(a-{t n )a + {t' n )..a-(t 1 )a + {t' 1 )), (B2) 


where p(0) = \g) (#|® |0) b ( 0 | and C is the Lindblad super¬ 
operator of the system, defined in the main text. C con¬ 
tains a deterministic part, which generates an evolution 
driven by H e ff and which conserves the number of exci¬ 
tations, and a jump part, which reduces the number of 
atomic excitations. Because of the form of the correla¬ 
tors, which contain an equal number of atomic creation 
and annihilation operators, the jump part of the evolu¬ 
tion of the operators gives a vanishing contribution to 
the correlation function, proving what was stated above. 


which, expressing explicitly the time dependence, is equal 
to 


(cr_e lHeff ( tn ^V + ..cr_e lHe{{( ^ tl t iV + ). 


(B3) 


Since H e ff is diagonal in the | g) , \e) basis, we can insert 
identity operators between the a operators in the form of 
| g) (g | + |e) (e\ in order to evaluate the Hamiltonians in the 
exponents. We immediately end up with n7=i 
where we ha ve d efined a = uj eg — iT/2. Inserting this 
result in Eq. (Bl) we get 


Appendix B: Examples of S-matrix elements 

In this section of the Appendix we provide explicit ex¬ 
amples of how to use the results derived above to calcu¬ 
late specific matrix elements. In particular we present 
the case of the scattering of 1 ) n photons on a two-level 
atom coupled to a one-directional waveguide and 2 ) one 
and two photons scattering on a Rydberg-EIT system. 


^( n ) _ 


(-0 


n /* + oo pl n pt 1 

/ dt n / dt ' n ... / dt[ 
J— oo J — oo J — oo 


[ p;k l (27r) 

n 

JJ g i ( Pi -a)ti e -i(ki-OL)ti 


+ permutations 




(B4) 


where the permutations are over the two sets of incoming 
and outgoing photon frequencies ki and pi . The integral 
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gives 


*^~K£)'{n (s>r* 

n / m— 1 \ —1 n 

X n km CM “1“ E A i J + permutations > x 

m= 1 ' 2=1 ' ^ 

n 

x ^(E A*)i (B5) 
2=1 

where = ki — pi, which coincides with the result of 
Ref. [22]. 


2. One and two-photon scattering on a 
Rydberg-EIT system 


The effective Hamiltonian of the Rydberg-EIT system 
is given by 

#eff = #0 + #hc + 7^ 5-] (B6) 

where i^o is given by Eq. ( [IT] ), and the hardcore interac¬ 
tion is 

Hue = y E( a ] a i + s ] s j)( a ] a j + s ) s i ~ !) ( B7 ) 

j 


with Uo oc corresponding to the three-level atom. 

For the single incident right moving photon with mo¬ 
mentum k , it follows from Eqs. © and © that the S- 
matrix element is 


with 


c(l) 



dtdt j (pt-kt') 
7T 


x(o|ro+(t)ot(i')|o), 



(B8) 


(B9) 


This element describes the amplitude of a single trans¬ 
mitted photon with momentum p, using the time ordered 
correlation function of the system operator. We set for 
notational simplicity the speed of light c = 1. The sec¬ 
ond term of S^-k+t he., the Fourier transform of the 
correlator 

/ -Too jjj .jp 

— e i(pt ~kt‘) ( 0 | Td+(t)dUt') | 0 ) 

-oo ^ 

= E e_i(fein+fc)(Zi ' Zj) [ G o(fc)]“; (BIO) 

b' 


can be obtained from the Green’s function Go(uj) = 
1/(cj — Hi) with the single-particle Hamiltonian 

H , = ( (-v - i£)<*« - IMy \ 

1 V 0 / 

(BH) 

Here, we have expressed Hi in the basis {|a^), |s^)}, 
and [G^k)]?? denotes the element (cq|Go(&) |cr') with 

cr, a' = a, s. The element = T^S p k gives rise to 

the transmission coefficient 

T k = 1 - i T -y Y}G 0 {k)]i*e- ikUzi ~ Zi) - (B12) 

*3 


For two right-going incident photons with mo¬ 
menta ki and k<^ the two-photon ^-matrix element 
(2 s ) 

P 2 j t -k 1 + k 2 + describes the amplitude to a final state 
with two transmitted photons of momenta pi and P 2 . By 
Eqs. © and Q in the main text, we find that 


*^Pl + ,P2 + ;/Cl + ,/c 2 + _ Tk 1 Tk 2 (Sp 1 k 1 S p 2 k 2 + ^Pl/C2^P2fel)(B13) 


r? 


ID 


EE e ifein 


4(2tt) 2 


“T-2^2 ) 


g: 




(pi,p 2 -,k 1 ,k 2 )e ikia(z n+ z i2) ^ 


with 


G aa:aa / 7 7 \ 

ii4;jU2(Pl’P2;fci,fc 2 ) 


= J dt[dt , 2 dt 1 dt 2 e^ Plt i +p <- kltl - k2t2 ' 1 

x (Tfii! (ii)a i2 (4)al Ti)aL (* 2 ) 


(B14) 


Here, ( Ta,, (t) )d,i 2 (t' 2 )aj i Ti )®J 2 fep denotes the con- 

nected four-point Green’s function, which involves 
an interaction between the two excitations in the 
system. The S-matrix also contains terms aris¬ 
ing from disconnected correlation functions, e.g., 

(VdqfTi)^^!)^> ^ Tdi 2 (t' 2 )a} j2 (t 2 describing the linear 
propagation of ea ch excitation separately, which yields 
the term in ( B14| ) proportional to T^T^. 

The interaction between two excitations under the 
Hamiltonian -Hhc+X^ Uijs\sjSjSi/2 can be represented 
by the ladder diagram shown in Fig. 6. This diagram is 
represented mathematically by the two-body T-matrix, 
T(E), which satisfies the Lippmann-Schwinger equation, 




T(E) = U + XJU 0 (E)T(E). 


(B15) 


Here, E = k\ + k 2 is the total energy, the vacuum bubble 
n 0 (E) = (E — H 2)- 1 is given in terms of H 2 = Hi (g) 
I 2 n + Hn 0 Hi, and the interaction matrix U has the 
diagonal element U c - f = = U-- 1 = Uo and Uff = Uij 

in the basis {| a^aj ), | diSj ), | s^aj ), The solution 

of Eq. ( |B15| is T (E) = l/fTT 1 ~n 0 {E)). The S-matrix 
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FIG. 5. The Feynman diagram for the T-matrix. 


We compare the result of g^ 2 \x) from the scattering 
theory and that from numerically solving the effective 
spin model (12) with the weak driving field in Fig. [4^ in 
the main text, which shows that they agree with each 
other perfectly. 


Appendix C: Generalized master equation 

In this section, we extend our derivation of the atomic 
master equation, in order to calculate the response of 
the system to an arbitrary few-photon input state (as 
opposed to a classical or coherent state). The initial state 
is generally written as |^in) = |<An) 0 |Xin), where |xin) 
is the initial state of the system, and 


can then be written as 

( 2 ) 

*^Pl + ,P2 + ;fcl+,fc 2 + = Tk 1 Tk 2 5p 1 k 1 fip 2 k2 — l -^-bpi+P2M+k 2 x 


*7* 3 

<7 1 <7 1 ]<j 2 cr 2 

+ (pi P2 )f 

where 


(B16) 


Mki,k 2 )]y; = [Gofa)]?,?. 

3132 

(B17) 

The Fourier transform of + . ki+ k2+ results in 

the wavefunction 


|^in) — ^ 1 (Cl) 

{n v ,k} k 

describes an arbitrary state of incident photons in the 
Fock state basis by the wavefunction ^in({^/c})- By the 
relation 




(C2) 


between the Fock state | n) and the coherent state | J) = 
J2n Jn I n)/Vn\, we rewrite | ip- m ) = F\ {J v ,k}) by the op- 
erator 




lim 

{Jv,k }~^0 


El ( / 9 in({ n v,/c}) 

{ n k } k 


1 £) n v,k 

\AFfe! 9JvT 


(C3) 


Jj( x x) — [ dpidp2 o(2) i(piX! 

ip{X c ,X) J ^ ^ > PiP2',k 1 k 2 e 


+P2X 2 ) 


0 \Ex c 


27r 


-{ 2T l+fc T f-fe' 


<kx)~ i-f J2 


aia 1 

*j 


*3* 3 

<j i<y \ \ < j 2 <j 2 


x[T(£)]^r CT >(§ + k, f- k%f}, 


(B18) 


of two transmitted photons, where the relative momen¬ 
tum k = (ki — & 2 )/ 2 , the center of mass coordinate 
x c = (xi+X 2 )/ 2 , and the relative coordinate x = X\—X 2 - 
The symmetric function 


, „ e -i k ixi { Zil +z i2 ) 

l F ( x W =-' 1 ^2-Er———^2xi(a il )xi'{a i2 ) 

1 V ll' 


xU^)xH^)[e^ 

+(x —x) 


*1*2 
(f -£v)x 


6{x) + e 


-i(^| -£i)x 


9{-x)] 


(B19) 


is defined by the eigenstates | xi) and | xi) of T~Li and H\ 

with the corresponding eigenenergies Si and e z *, (cq |%/) = 

a nd (o -j \xi) = Knowledge of the wave 

function ( B18| ) in the case where the two incident photons 

have the same momentum, k\ = k 2 = E/2, enables one 

to calculate the second-order correlation function for the 

2 

outgoing field, g^ 2 \x) = tt^(x c ,x ). 


acting on the coherent state \{Jk})- 

The evolution of the reduced density matrix p s (t) = 
Trb[U(t)p(0)W(t)\ is determined by 

U(t) = Tex p[—i f dsH(s )], 

Jo 

where p(0) = IV’in) (V’inl is the density matrix of the initial 
state, and H(t) is the total Hamiltonian in the rotating 
frame with b Vj k —^ b v ^e~ lvkt . By the relation \(p[ n } = 
J 7 \{Jv,k}}, the reduced density matrix reads 

Ps(t) = J-*J-Tr b [U(t) I {J v , k }) ({J v , k }\ Ps(0)U\t)}, (C4) 

where p s (0) = |xin) (Xin| is the initial density ma¬ 
trix of the system. The displacement transformation 
Vd \{J k }) = e2^ v ’ k \ Jv ’ k \ |0> leads to 

P.(t) = P*Pe^P J ^' 2 pj(t), (C5) 

and the generating density matrix 

pj{t) = Tr b [U d {t) |0> (0| p s ( O)Wl(t)], (C6) 

where the unitary transformation 

U d (t) = T exp[-i [ dsH d {s )] (C7) 

Jo 
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1 ph inside system: 

No EIT resonance 

-* 

2 ph inside system: 
EIT resonance 

2 ph outside system 

f 

1 ph inside system: 
No EIT resonance 


o l y 


FIG. 6. An intuitive explanation of how the two-photon wave 
function evolves within the atomic medium (located in the re¬ 
gion 0 < x,y < L) can be found by considering a larger space, 
where first the two photons are both outside the medium 
(x,y < 0) and where one photon enters the medium first. 
We specifically consider the case where infinite-range inter¬ 
actions cause the two photons to align with the EIT trans¬ 
parency condition when they are both inside the medium, as 
discussed further in the main text. The two-photon trans¬ 
mission is proportional to the value of the two-photon wave 
function at x,y = L. 

is given by H^(s) = VdH(s)V^. The Hamitonian 
Ha(s) = H(s ) +Hj(s ) is obtained by replacing the oper¬ 
ator b V: k by b v ^ + J v ,k in the Hamiltonian H(s), where 
Hj(s ) describes the system under the driving fields J V: k- 
In the density matrix pj(t), the initial state becomes the 
vacuum state due to the displacement transformation, 
thus the evolution of pj(t) satisfies the master equation 

d t pj(t) = Cpj(t) - i [Hj(t),pj(t)\, (C8) 

where C is the Lindblad super-operator of the system 
without Hj. 

In conclusion, Eq. ( |C5| ) establishes the relation bew- 
teen the evolution of the reduced density matrix p s (t) for 
a few incident photons (which could be a non-classical 
state) and the evolution of the reduced density matrix 
pj(t) for the system with classical (coherent state) driv¬ 
ing fields J v ,k- As a result, the few photon scattering 
problem can be understood as the perturbation expan¬ 
sion of the driving strength J Vj k- 

Appendix D: Linear approximation for two photon 
transmission T 2 

In this section of the Appendix we provide an intuitive 
explanation based on linear optics to the behavior of the 
two-photon transmission T 2 depicted in Fig. UK 13 - % 
indicating with x and y the coordinates of the first and 


second photons, we can divide the space into four regions 
according to their positions: both photons outside the 
atomic medium, one photon inside and one outside, and 
both photons inside the medium, see Fig. [6j As we are 
studying an infinite-range interaction, within each region 
the evolution is effectively linear. The different disper¬ 
sion relations in each region, however, leads to non-trivial 
boundary conditions at their edges. Relevant to our dis¬ 
cussion is the value of the two-photon wave function at 
the boundaries 0 < x < L and 0 < y < L. In particular, 
evolution in the region 0 < x,y < L with these boundary 
conditions dictates the two-photon transmission, which 
is proportional to the value of the two-photon wave func¬ 
tion at (x,y) = (L, L). 

We study the case where the two photons have detun¬ 
ings A = C/2. In this case, infinite-range interactions 
cause these two photons to satisfy the EIT transparency 
condition when both photons are inside the medium. 
Then, there are two qualitatively different regimes for 
the boundary values of the two-photon wave function, 
depending on the detuning from two-photon resonance S 
(for simplicity we assume Sl = 0, so that A = 5—Sl = S): 
the small detuning regime, S < A eit, in which a single 
photon can travel with high transmission through the 
medium and reach the detection point; and the large de¬ 
tuning regime, 5 > A eit, in which a single photon is 
absorbed and cannot reach the detection point. 

Small detuning regime. — Within the first regime S = 
C/2 < A eit, a single polariton still fits within the 
EIT transparency window and exhibits high transmission 
through the medium. A representative probability distri¬ 
bution for a single excitation in this regime is illustrated 
in Fig. [7^. Now, we consider the two-excitation manifold. 
After detection of the first outgoing photon, the second 
photon can be at any position inside the ensemble with 
nearly uniform probability, and on average has to travel 
over a distance Nd/2 before reaching the detection point. 
Therefore we expect T 2 = | exp(— ik(S)Nd/2)\ 2 = \fT[. 
The predicted behaviour is in good agreement with full 
simulations results as shown in Fig. [8] for small detuning. 

Large detuning regime. — In this regime a single photon 
is strongly absorbed. This is illustrated in Fig. [ 7 ) 3 , where 
one observes a strong decay of the single-polariton prob¬ 
ability and field intensity inside the medium. Within the 
context of Fig. [6j this localized single excitation serves 
as the boundary condition on the segments 0 < x < L 
and 0 < y < L. Furthermore, the evolution in the re¬ 
gion 0 < x,y < L is effectively linear. Thus, the effect 
of this boundary condition at the detection point (a two- 
particle problem) can be mapped onto a simpler problem, 
wherein one studies how a single excitation, initialized in 
the shape of the localized photon, transmits through the 
medium. In this regime, we therefore can estimate that 
T 2 is twice the value of the maximum transmission asso¬ 
ciated with this single localized photon, whose evolution 
we calculate numerically. 

The good agreement of our simplified estimates with 
full simulations is demonstrated in Fig. [8j in which T 2 as 
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(b) 



j 



C 

FIG. 8. T 2 versus C — 26: comparison of the full simulation 
T 2 (red triangles) with the simple models for the small de¬ 
tuning regime (where one expects T 2 = \/T\, black line) and 
the large detuning regime (blue curve), where the two-photon 
transmission is obtained from the propagation of a single pho¬ 
ton localized at the front of the medium (near z = 0). 


FIG. 7. Single photon probability (bo Ut (jd)b ou t(jd)) (green) 
and single polariton probabilities {cr° ee ) (black) and ( a J ss ) (red) 
as function of atom position j, for asymptotically large time 
at A = 6 = 0.1 (a) and at A = 8 = 1 (b). Probabilities are in 
unit of a 2 . Other parameters: N = 80,11 = 2, F = 3, Fid — 
l,cx = 10 -6 . 


a function of C = 26 is compared for the different cases 
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